N90- 25055 


Minimizing Distortion and Internal Forces in Truss Structures by 

Simulated Annealing 

by 

Professor Rex K. Kincaid 
Mathematics Department 
The College of William and Mary 
Williamsburg, VA 23185 


Inacuracies in the length of members and the diameters of joints of large truss reflector 
backup structures may produce unacceptable levels of surface distortion and member forces. 

However, if the member lengths and joint diameters can be measured accurately it is possible to 
configure the members and joints so that root-mean-square (rms) surface error and/or rms 
member forces is minimized. 

Following Greene and Haftka (1989) we assume that the force vector f is linearly 

proportional to the member length errors of dimension NMEMB (the number of members) and 

joint errors e of dimension NJOINT (the number of joints), and that the best-fit displacement 

vector d is a linear function of f. Let NNODES denote the number of positions on the surface of 

the truss where error influences are measured. Let U. , (NNODES x NMEMB) and U r (NNODES 

M J 

x NJOINT) denote the matrices of influence coefficients. Then d = U. .e. _+ U T e T . Concatenating e x 

M M J J M 

with e and with U f yields d = Ue. 

Let D be a positive semidefinite weighting matrix (in our computational experiments we let 
D be an identity matrix) denoting the relative importance of the surface nodes where distortion is 
measured. The mean-squared displacement error can then be written as 


2 T T T 

d = e U DUe = e He. 
rms 


£ 

A similar construction can be derived foj mean-squared member force error, s (see Greene and 

Haftka (1989)). Minimizing d (or s ) can be formulated as a combinatorial optimization 
, , _ _ , . , r ms . r ms. . . , , . . .2 

problem. That ts, finding the permutation of the components of e and e that minimizes d (or 
2 . , : ... .2 . 2 M ,J , rms 

s ) is equivalent to minimizing d (or s ) directly. Unfortunately there 

(W?JfEMB!)*(NJOINT!) possiblities r to consider. 51 Hence, an enumeration scheme is out of the 

question. However there are many combinatorial optimization problems with exponentially large 

solution spaces that can lie solved by algorithms whose time complexity is bounded by a 

polynomial function of die problem parameters. 

To classify this problem we compare it to a similar combinatorial optimization problem. In 

particular, when only the member length errors are considered, minimizing d is equivalent to 


rms 
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the quadratic assignment problem. The quadratic assignment problem is a well known NP- 
complete problem in the operations research literature. Hence, minimizing df is also an NP- 
complete problem. Moreover, if a problem is NP-complete it is highly unlikely that an algorithm 
exists which can determine an optimal solution in polynomial time and, therefore, (polynomial 
time) heuristic solution techniques should be employed. Greene and Haftka (1989) tested two 
heuristics of the same type. Thpy use pairwise interchange and triple interchange of the 
members and joints to reduce d f . The fjpeus of our research has been the development of a 
simulated annealing algorithm to reduce d . The plausibilty of this technique has been its 
recent success on a variety of NP-complete combinatorial optimization problems including the 
quadratic assignment problem. 

Simulated annealing was first proposed and used in statistical mechanics in the early 
1950’s (see Metropolis et al. (1953)). However, not until Cemy (1982) was simulated annealing 
used to solve a NP-complete combinatorial optimization problem— the traveling salesman 
problem. A physical analogy for simulated annealing is the way liquids freeze and crystallize. 

As the liquid is cooled slowly the atoms line themselves up and form a pure crystal that is 
completed ordered. The pure crystal is the minimum energy for this system. The basic 
procedure consists of a loop over a random displacement generator that produces changes in the 
objective function value. If this change is negative the displacement is accepted and the 
objective funedon is reduced. If this change is non-negative the displacement is accepted 
probabalistically. That is, uphill climbs are accepted with some positive probability which 
decreases as the temperature decreases. Simulated annealing must be used with some care. In 
addition to determining how to generate random displacements, one must also pick a starting 
temperature T, a cooling rate TFACTR, and a stopping temperature T^. If these parameters are 
not chosen appropriately simulated annealing may produce poor results and/or run for an 
exponential amount of time. 

Figure 1 is a graph of the objective function value (d ) for ten random starting 
arrangements of the components of e for three different heuristics. All computational 
experiments were done on a MicroVAX. The two interchange heuristic is very fast (an average 
epu time of 1 . 1 minutes per run) but produces widely varying results. The two and three 
interchange heuristic provides less variability in the final objective function values but runs much 
more slowly (an average epu time of 68 minutes per run). Simulated annealing produced the best 
objective function values for every starting configuration and was faster than the two and three 
interchange heuristic (an average epu time of 42 minutes per run). 
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C ********** MAIN PROGRAM ************* 

C 

INTEGER NMEMB, N JOINT, NROW 

PARAMETER (NMEMB-102 f N JOINT-3 1, NROW^NMEMB+N JOINT NW-19) 
INTEGER IORDER (NROW) , NDIM ' 

DOUBLE PRECISION IT, M, ZQAP, T f TFACTR, ZCHX, PSUM, TSUM 
DIMENSION H (NROW, NROW) , M<NROW) , PSUM (NROW) 

$, DFDE {NW, NROW) , DWDE (NW, NROW) 

CHARACTERS 5 MSG 
REAL TIM (20) 

OPEN (UNITES, FILE-' QAP2 . IN' , STATUS-' OLD' ) 

OPEN (UNIT-6, FILE-' QAPJNT.OUT' , STATUS-' UNKNOWN' ) 

OPEN (UNIT— 7, FILE—' MTEN . DAT' , STATUS—' OLD' ) 
c 

^*******^*************^^****^****^^*^ 

c 

C Raid in input data. Influence matrix H-UDU, member 

C length error* M, joint diameter error* M, displacement 

C derivatives DWDE, force derivatives DFDE, and initial 

C objective function value ZQAP. The input file QAP . IN 

C is created by GENQAP . FOR . 

C 

C «********0*0*MO***00**0**t*M**0#0*«*****#*** # * ###t * 

C 

DO 21 I— 1, NROW 

READ(5, 901) (H (I, J) , J— 1, NROW) 

21 CONTINUE 

DO 20 1-1, NMEMB 

READ (5, 901) (DFDE (I, J) , J— 1, NROW) 

20 CONTINUE 

DO 22 1*1, NW 

READ (5, 901) (DWDE (I, J) , J-l, NROW) 

22 CONTINUE 

DO 2400 J-l, 3 

DO 17 I— 1 , NROW 

IORDER (I) -I 

17 CONTINUE 

READ (7, 901) (M ( I) , I— 1, NROW) 


READ ( 7 , 902 ) ZQAP 
C 


READ (7, 900) MSG 


C 

C** 

C 

c 

c 

c 

c 

c 

c** 

c 


Use the largest eigenvalue of H to provide a bound on 
the difference between the largest and smallest objective 
function v»lu«». For thi» H, 9.779335 i« th« 4pproprl«t« 
eigenvalue. r 
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900 

901 

902 


T-0 .0 

DO 79 I— 1 , NROW 

T— T+M (I) *M(I) 

CONTINUE 
T— T*10*9 .779335 
TFACTR-0 . 96 
FORMAT (IX, A) 

FORMAT (IX, 5E16. 12) 
FORMAT (IX, E16. 12) 


c ********* End initialization and echo results 

a 

a 

WRITE (6,*) 'ITERATION ',J 

WRITE (6,*) 'Start Temperature- T, TFACTR 

WRITE (6,*) 'Starting ZQAP-' , ZQAP 

CALL SECOND (TIM (J) ) 

C 

CALL ANNEAL (M, H, IORDER, NMEMB, NROW, TFACTR, ZQAP, T) 
CALL SECOND (TIM(10+J) ) 

WRITE ( 6, *) 'Execution time ' , TIM(10+J) -TIM(J) 

„ WRITE (6,*) 'Final annealing objective value ' , ZQAP 

CALL OBJCHK(M, IORDER, H, PSUM, NROW, ZCHK) 
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WRITE (6,*) 'Obj. value check. ' , ZCHK 
WRITE (6, *) (IORDER (I) , I-1,NR0W) 

WRITE (6,*) 'Final temperature ',T 

2400 CONTINUE 
STOP 
END 
C 
C 

c 

SUBROUTINE ANNE AL (M, H, IORDER, NMEMB, NROW, TFACTR, ZQAP, T) 

C 

c 

c ************************************ ********************************** 

c 

c Thia algorithm find* the permutation of tha componanta 

C of tha vactor M that minimi zee tha product MHM for any 

C raal ayimnatria poaitiva dafinita matrix H. Tbara ara 

C NROW ( NMEMB +N JO I NT) componanta of M and H ia NROW by NROW. 

C Tha array IORDER (I) apacifiaa tha pannutation of M. On 

C input, tha elamanta of IORDER may ba aat to any permutation 

C of tha numbara 1 to NROW. Thia routina will ratum tha baat 

C altarnativa permutation it can find. 

C 

c T ia tha currant temparture. 

C NOVER ia tha max number of awapa triad at any tamparature T. 

C NLIMIT ia tha max number of auccaaaful awapa before continuing. 

C TFACTR ia tha annealing achedule, Tnew-T old* TFACTR. 

C ZQAP denote a tha objective function value at any time T. 

C DE danotaa tha change in ZQAP whan two componanta ara awappad. 

C 

C***** ***************** * *********************************************** 

c 

c 

INTEGER NMEMB, IORDER (NROW) , N(2) , NOVER, NLIMIT, IDUM 
DOUBLE PRECISION M, FT, TFACTR, ZQAP , DE, T, TSUM 
DIMENSION M < NROW) , H (NROW, NROW) 

LOGICAL ANS 
NOVER* 10*NROW 
NLIMIT-l*NROW 
IDUM*- 1 
NSUCC*1 
NCNT*0 

N JO I N T -NROW - NMEK1 
C 

C******* Loop until tamparature ia too email or NSUCC-0 . 

C 

DO WHILE (NCNT.LT. 600 .AND .NSUCC.GT.O) 

NCNT-NCNT+1 

NSUCC-0 


C 

C******* Local aaarch of neighbor* of currant assignment 
C 

DO 12 K— 1, NOVER 
C 

C* ****** N(l) and N{2) ara tha two component* of M to ba awappad. 
C 

IF (RAN3 (IDUM) . GT . 0 .76692) THEN 

N (2) -1+INT (NJOINT*RAN3 (IDUM) ) 

N (1) -1+INT ( (NJOINT-1) *RAN3 (IDUM) ) 

IF (N(2) .EQ.N(l) .AND.N(2) .EQ.NJOINT) THEN 
N(1)*N(1)-1 

ELSE IF (N(2) .EQ.N(l) ) THEN 
N(2)-N<2)+1 

END IF 

ELSE 

M (2) -1+INT (NMEMR*RAN3 (IDUM) ) 

N ( 1) —1+INT ( (NMEMB-l) *RAN3 (IDUM) ) 

:CF (N(2) .EQ.N(l) .AND. N(2) .EQ. NMEMB) THEN 
N(l)-N(l)-1 

ELSE IF (N{2) .EQ.N(l) ) THEN 
N(2)-N(2)+l 

END IF 

END IF 

CALL SWPCST(M,H, IORDER, NMEMB, NROW, N, DE) 
CALL METROP (DE, T, ANS) 

IF (ANS) THEN 

NSUCC-NSUCC+1 

ZQAP -ZQAP +DB 

CALL SWAP ( IORDER, NROW, N) 
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IF (NSUCC.GE .NLIMIT) GOTO 2 


12 CONTINUE 

2 T»T*TFACTR 

END DO 

WRITE (6, *) ' NCNT' , NCNT 

RETURN 
END 
C 


SUBROUTINE SWPCST (M, H, IORDER, NMEMB, NROW, N, DE) 

C 
C 

C**************************************************** 

c 
c 
c 
c 
c 
c 

£♦* *********** ********************************** ********************** 

c 

c 


Thia aubrout in* return* tha valua of tha change in tha 
objactive function for a propoaad swap of two position* 
in tha currant parmutatlon aaaigrunant IORDER. On output 
DE ia tha valua of tha changa {+ or -) . 


INTEGER NMEMB, IORDER (NROW) , N (2) , II, J1,K,K1, ITMP 
DOUBLE PRECISION M, H f LTSUM, RTSUM, DIFF, DE, SQDIFF 
DIMENSION M (NROW) , R (NROW, NROW) 

C 

c ********* initialization 
C 


DE=0 .0 
RTSUM=*0 . 0 
LTSUM^O .0 
Il=IORDER (N (1) ) 

Jl" IORDER (N (2) ) 

C 

C********* p Ut indicaa of M in aacanding ordar, 
C 


II < J1 


IF (Il.GT.Jl) THEN 
ITMP-I1 
NTMP-N(l) 
Il-Jl 
N(1)*N<2> 
Jl-ITMP 
N(2)«NTMP 

ENDIF 

C 


c********************************** # * ##tk#jk###A##A### 


******************** 


c 

c 

c 

c 

c 

c 

c 

c 

c< 

c 


Thia aaction of tha coda computa* tha changa in tha objactiva 
function valua, DE, in linaar tima. To do thia, a pointar array 
IORDER la uaad to kaap track of tha awltchaa in tha array M. 
Sinca only two components of M ara awitchad at any ona tin* only 
two rowa and two columns of tha matrix H naad ba conaidarad to 
computa DE . 


DO 12 K*l,NROW 

Kl»IORDER (K) 

IF (Kl.EQ.Il.OR.Kl.EQ. Jl) GOTO 12 

LTSUM-LTSUM+H (K, N < 2) ) *M (K1 ) 
RTSUM-RTSUM+H(X,N<1) ) *M(K1) 

12 CONTINUE 

DIFF*=M ( Jl ) -M ( II) 

SQDIFF" (M< Jl) ** 2 ) - (M(I1) **2) 

DE»(SQDIFF*H (N(l) , N(l) ) ) + (2*DIFF*RTSUM) 

$ - (SQDIFF*H(N{2) ,N(2) ) ) - (2*DIFF*LTSUM) 

RETURN 
END 
C 

c 


c 

c 

e 

c 

c 

c 

c 

c 


SUBROUTINE SWAP (IORDER, NROW, N) 


Thia routina parforma tha actual awap in IORDER batwaan 
poaitiona N(l) and N(2). On output IORDER ia modifiad to 
r*flact thia axchanga. 


C******* *********** ************************** 

c 
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c 

INTEGER NROW, I ORDER (NROW) , N (2) , ITMP 
C 
C 

ITMP-IORDER<N(l) ) 

IORDBR(N(l) )*IORDER(N(2)) 

IORDER(N(2) )-ITMP 
RETURN 
END 
C 
C 

S UBROUT INE METROP ( DE , T , ANS ) 

C 

C 

C** ******************* ** ******* **************** ******* ********** 

c 

C Matropolia algorithm. ANS 1 a a logical variable which 

C iaauaa a vardict on whathar to accapt a raconf iguration 

C which laada to a changa DE in tha objactiva function B. 

C If DE<0 , ANS - .TRUE., whila if DE > 0, ANS ia only 

C .TRUE, with probability axp(-DE/T), whara T ia a 

C tamparatura datarminad by tha annaaling achadula. 

C 

c **** ******************* **************************************** 

c 

c 

DOUBLE PRECISION DE, T 
PARAMETER ( JDUM-1 ) 

LOGICAL ANS 

ANS* (DE . LT .0.0) .OR. (RAN3 ( JDUM) .LT .EXP (-DE/T) ) 

RETURN 

END 

CC 

C 

FUNCTION RAN3 ( IDUM) 

C 

C 

C* ************************************************** ************** 

C 

C Returns a uniform random davlata batwaan 0.0 and 1.0. 

C Sat IDUM to any nagativa valua to initlallza or 

C raintializa tha aaquanca. (see Numarical Racipaa p. 199) 

C 

c ****** *********************************************** ************ 

c 

C 

PARAMETER (MBIG-1000000000, MSEED-161903398 , MZ-0, FAC-1 ./MBIG) 
DIMENSION MA (55) 

DATA IFF /0/ 

'*****InitialiTat Lon 

IF (IDUM.LT. 0 .OR. IFF.EQ.0) THEN 
IFF-1 

MJ-MSEED- IABS (IDUM) 

MJ-MOD (M.T, MBIG) 

MA (55 ) “Mi/ 

MK-1 

DO 11 1-1 , 54 

n-MOD(21*I,55) 

MA { I I ) -MX 
MK-MJ-MX 

IF (MK . LT . M2 ) MX-MK+MBIC 
MJ-MA(II) 

11 CONTINUE 

DO 13 K-1,4 

DO 12 1-1,55 

MA (I) -MA (I) -MA(1+M0D (1+30, 55) ) 

IF (MA (I) . LT.MZ) MA(I) — MA{I) 4-MBIG 

12 CONTINUE 

13 CONTINUE 
INE XT-0 
INEXTP-31 
IDUM— 1 

END IF 

C 

c ********* End initialization 
C 

INEXT-INEXT4-1 

IF ( INEXT . EQ .56) INE XT- 1 
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C 

C 


c 

c 

c 

c 


c 


4 


5 


C 

c 


INEXTP- INEXTP + 1 

IF ( INEXTP . EQ .56) INEXTP-1 

MJ*MA ( INK XT ) -MA ( INEXTP ) 

IF (MJ . LT .M2) MJ-MJ+MBIG 

HA (IKE XT) =MJ 

RAN3*MJ*FAC 

RETURN 

END 


SUBROUTINE SECOND (TIM) 
TIME0=*0 . 0E+00 
TIM^SECNDS (TIMEO) 
RETURN 
END 


SUBROUTINE OBJCHK (M, I ORDER, H, PSUM, NROW, ZCHK) 

INTEGER NROW, I ORDER, II, J1 
DOUBLE PRECISION ZCHK, M, H, PSUM 

DIMENSION M (NROW) , H (NROW, NROW) , I ORDER (NROW) , PSUM (NROW) 

ZCHK=»0 .0 
DO 5 1*1, NROW 

Il**IORDER ( I ) 

PSUM ( I ) *0 . 0 
DO 4 J»l, NROW 

Jl-IORDER (J) 

PSUM (I) «=PSUM (I) +H (I, J) *M ( Jl) 

CONTINUE 

ZCHK - ZCHK + PSUM ( I ) ♦M (II) 

CONTINUE 

RETURN 

END 
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